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ABSTRACT 

We analyse the dynamical properties of substructures in a high-resolution dark matter simu- 
lation of the formation of a Milky Way-like halo in a ACDM cosmology. Our goal is to shed 
light on the dynamical peculiarities of the Milky Way satellites. Our simulations show that 
about 1/3 of the subhalos have been accreted in groups. We quantify this clustering by mea- 
suring the alignment of the angular momentum of subhalos in a group. We find that this signal 
is visible even for objects accreted up to z ~ 1, i.e. 8 Gyr ago, and long after the spatial coher- 
ence of the groups has been lost due the host tidal field. This group infall may well explain the 
ghostly streams proposed by Lynden-Bell & Lynden-Bell to orbit the Milky Way. Our analy- 
ses also show that if most satellites originate in a few groups, the disk-like distribution of the 
Milky Way satellites would be almost inevitable. This non-random assignment of satellites to 
subhalos implies an environmental dependence on whether these low-mass objects are able 
to form stars, possibly related to the nature of reionization in the early Universe. With this 
picture, both the "ghostly streams" and the "disk-like configuration" are manifestations of the 
same phenomenon: the hierarchical growth of structure down to the smallest scales. 

Key words: methods: numerical - Galaxy: formation - galaxies: dwarf - galaxies: kinematics 
and dynamics - dark matter. 



1 INTRODUCTION 

In cold dark-matter (CDM) cosmologies, large galaxies are the re- 
sult of the aggregation of smaller subunits. Some of these subunits 
may survive until the present day in the form of satellites, while 
some may be completely destroyed in the course of time, and con- 
tribute to the field. In A^-body CDM simulations of the formation 
of galaxy-size dark-matter halos, the satellites (often referred to 
as substructures or subhalos) show significantly different proper- 
ties than those of the "luminous" satellites around galaxies like the 
Milky Way. An example of this discrepancy is the "missing satellite 
problem": subhalos outnumber the bright satellites by factors of a 
hundred or more ( Kauffmann, White & Guiderdoni|199"3| |Klypin| 
|et al.||1999[ |Moore et al. 1999). Furthermore, their spatial distri- 
bution is typically much shallower than observed for the luminous 
satellites jGao et al.|20"04a||Taylor, Babul & Silk|2004^ . Therefore, 
the relation between resolved substructures/subhalos in dark mat- 
ter simulations and the luminous satellites in galaxy halos is still 
unclear. Attempts to reconcile these two populations, using semi- 



analytic models of galaxy formation ( Kauffmarm, White & Guider-| 


|doni|1993||Benson et al.|2002| 


Kravtsov, Gnedin & Klypin|2004 1 or 


full-fleshed SPH-simulations 


Maccio et al.|2006||Libeskind et al.' 



|2007 ) have produced interesting results, and helped us gain insight 
into the relevant processes on the smallest galactic scales. 



In the past ten years new attention has been drawn to the 
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properties of the satellite population in the Local Group. Starting 
with Lynden-Bell & Lynden-Bell' ('1995'), the existence of ghostly 
streams of satellites (dwarf galaxies and globular clusters) was pro- 
posed. These objects would share similar energies and angular mo- 
menta producing a strong alignment along great circles on the sky 
(Palma, Majews ki & Iohnston|2002"| . Recently, it has been argued 
that the MW satellites define a disk-like structure, that is so highly 
flattened {rms thickness of 10 — 30 kpc) that may at first sight be in- 
consistent wifli CDM models jKroupa, Theis & Boily|2005||MetzIl 
|Kroupa & Jerjen|2007^ . However, sophisticated modelling combin- 
ing semi-analytic galaxy formation recipes with dark matter simu- 
lations has produced results that are consistent with observations 
l |Kang et al.|2005||Libeskind et al.|2005||Zentner et al.|2005| l. 

In the past two years, about a dozen low surface brightness 
dwarf satellite galaxies around the MW have been found in the 
Sloan Digital Sky Survey (|Willman et al.||2005aP [Belokurov et| 
|al.|2006|[2007l|Zucker et al.|2006a|b[|Irwiii et al.|2007> . This is an 
increase by a factor of two in the number of known satellites. Be- 
cause of the strong selection bias (SDSS focus is on the north galac- 
tic pole), it is unclear whether these satellites confirm the flattened 
disk-like structure. It is also unclear whether they help solve the 
missing satellite problem ( |Simon & Geha|2007| ). Satellites around 
M31 also show peculiarities in their distribution, such as an excess 
of objects on the side closest to the Milky Way ( [McCormachie & 
Irwin 2006 ), and possibly a similar degree of alignment ( |Koch &| 
|Grebelt2006^ 

These facts have motivated us to revisit the distribution and 
properties of the subhalo population in CDM simulations. In par- 
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Table 1. Numerical parameters for the GAnew-series simulations. The 
number of low-resolution particles remains at a roughly constant level of 
nlr ~ 1.272 X lO''. 



Name 


mp [/i-^Mo] 


Nhr 


Zstart 


e[h-^ kpc] 


GAOnew 


1.677x10* 


68 323 


70 


1.4 


GAlnew 


1.796x10^ 


637 966 


80 


0.8 


GA2new 


1.925x10*' 


5 953 033 


90 


0.38 


GA3new 


2.063x10^ 


55 564 205 


60 


0.18 



ticular, here we focus on the infall of substructures onto a Milky- 
Way like dark matter halo in a ACDM cosmogony. We use a high- 
resolution dark-matter simulation that is a variant of the GA series 
( |Stoehr et al.|2002[ |Stoehr 20061 as described in Section [TT] We 
find evidence of group infall onto the MW-like halo, which may ex- 
plain the ghostly streams proposed by| Lynden-Bell & Lynden-Bell| 
( |1995 1. This is discussed in Section |2.2| In Section |2.3| we show that 
the disky configuration of satellites is consistent with CDM if most 
satellites have their origin in a few groups. Thus both the "ghostly 
streams" and the "planar configuration" are manifestations of the 
same phenomenon: the hierarchical growth of structure down to the 
smallest galactic scales. We present our conclusions in Section|3] 



2 SUBHALOS IN THE DARK MATTER SIMULATIONS 
2.1 Description of tlie simulations 

We have analysed the GAnew series of high resolution simulations 
of a MW-like halo ( |Stoehr|2006| l. The simulations were carried out 
with GADGET-2 (Springel et al.'200 1J. The halo was sele cted from 
the M3 ACDM series (see Ciardi, Stoehr & White|2003^ with cos- 
mological parameters fio = 0.3, Q.a = 0.7, h = 0.7, and Hubble 
constant Ho=100/i km s~^ Mpc~^. Candidate halos in this sim- 
ulation were selected according to the following criteria: i) well- 
resolved, with > 500 particles; ii) environment similar to the Local 
Group (free from nearby rich galaxy clusters); Hi) peak circular ve- 
locity approximately 220 km s~^, and iv) no major merger since 
z = 2 (see |Gilmore, Wyse & NorTis||2002] >. After choosing the 
best MW-like halo, the zoomed initial conditions technique ( |Tor-| 
men, Bouchet& White|1997| l was used to produce a higher resolu- 
tion simulation. The MW-like halo was simulated four successive 
times with the mass resolution increased by a factor 9.33 each time. 
In the highest mass resolution simulation (GASnew), there are ap- 
proximately 10^ particles within the virial radius. Each of these re- 
simulations produced 60 outputs equally spaced in log (1/(1 + z)) 
between z — 37.6 and 2 = 0. Table [T] summarises the parameters 
and properties of the GAnew series. 

Virialised structures were identified in the high-resolution re- 
gion of the simulation using the standard friends-of-friends algo- 
rithm linking particles separated by less than 0.2 of the mean inter- 
particle separation. The SUBFIND algorithm |Springel et al.|200T) 
was then applied to each FOF group to find the gravitationally 
self-bound subhalos. The smallest resolved subhalo contains 10 
particles, and its mass is ~ 2.7 x 10^ Mq in GA2new and ~ 
2.9 X 10® M© in GA3new. The number of subhalos present in the fi- 
nal output is 504 for GA2new and increases to 3, 892 for GASnew. 
The virial mass of the MW-like halo is 2.4 x lO^^/i"^ M© and the 
virial radius is Rvir = 217 kpc at present time in GA3new. 



2.2 Group Infall of Dark Matter Substructures 

2.2.1 Infall pattern 

The bottom panels of Fig. [T] show the present-day spatial distri- 
bution of subhalos accreted at four different redshifts (zacc) in 
GA3new. We use the most bound particle of a subhalo to represent 
its position and velocity. The accretion epoch of a subhalo is de- 
termined as the time when its most-bound particle becomes part of 
the FOF group in which the MW-like halo is located. The reference 
frame used in Fig.[T]is defined by the principal axes of the MW-like 
halo. Its principal axes have been determined by diagonalising the 
inertia tensor 



(1) 



where xf is the i coordinate of the fith particle with respect to a ref- 
erence frame centred on the main halo, — (y^ )^ + (yi^/s )^ + 



iVs /'})''' ™d Vi its coordinates in the principal axes frame i Du 
[binski & Carlberg||1991| [Zentner et ari|2005| l. We only consider 
particles with ( ^ Rvir, and iterate until s and q have changed 
by less than ~ 10~^. The present time minor- to-major axis ratio of 
the MW-like halo is g = 0.6, and it has changed by less than 5% 
over the last 3 Gyr. 

The top row of Fig.[T|shows the spatial distribution of subhalos 
at the time they were accreted onto the MW-like halo. The (grey) 
dashed circles in top panels show the virial radius at that epoch. 
There are 789, 231, 133 and 117 subhalos which have survived un- 
til present time that were accreted at z = 0.06, 0.28, 0.63 and 1.08 
respectively. Note that in the first column, we plot for clarity only 
the first 200 most massive subhalos accreted at Zacc = 0.063. 

The distribution of subhalos is elongated roughly along the di- 
rection of the major axis of the MW-like halo at each epoch. The 
projected major axis orientations of subhalos accreted at a given 
epoch are indicated by the dashed lines in Fig. [T| These are com- 
puted using Eq. |T]( but without resampling The alignment be- 
tween the host and the subhalo's major axes is prominent except 
for the Zacc ~ 0.28 snapshot. At this time the accreted subhalos 
are much more isotropically distributed, preventing a clean deter- 
mination of the orientation of the principal axes of inertia (since the 
intermediate-to-major axis ratio is 0.98). 

Fig.[T]also shows that the distribution of subhalos is clumpy on 
small scales. These clumps are formed by subhalos sharing similar 
velocities (as shown by the arrows in this figure). This indicates 
that infall does not occur in isolation but in groups. This is very 
reminiscent of the way that clusters of galaxies grow through the 
mergers of groups, residing in the intersections of filaments (e.g. 
|Knebe et al.||2004l l. What we are seeing here is that also galaxy- 
size halos grow via accretion of "subgroup"-size structures. 

The clumps of subhalos seen in Fig. [T] share essentially the 
same angular momentum at the time of infall. This implies that, 
even if the spatial clustering is less prominent after a few orbits 
inside the main halo, the lumpiness may still be present in the space 
of angular momenta, provided these are nearly conserved. 

We quantify the degree of clustering by computing the two- 
point "angular correlation function", ui{a), of the present-day an- 
gular momentum of our subhalos. The angle a is defined by the 
relative orientation of the angular momenta of any two subhalos. 



^ When the number of subhalos is small, as in the right panel of Fig.^ an 
iterative procedure quickly leads to a sample with too few points, and hence 
to unreUable results. 
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Figure 1. Distribution of present-day sublialos in tlie Milky Way-like halo in GA3new. The top panels show their distribution at the redshift of accretion Zacc 
in the principal axis reference frame defined by the main halo at that epoch, while the bottom panels corresponds to the present-day. The arrows represent 
the velocity vector of each subhalo: its length coiTcsponds to 0.15 of the velocity magnitude while the orientation is defined by the direction of motion. The 
dashed, dash dot and dotted lines in each panel indicate the directions of the major, intermediate and minor axes of the ellipsoid defined by these subhalos. 
The circle denotes the virial radius of the host at the accretion epoch. Several groups of subhalos can be seen in these diagrams, in particular at the time of 
accretion. 



i.e. cosaij = • Lj/(|Li||Lj|). Therefore the correlation func- 
tion uj measures the number of pairs with otij < a compared to the 
expectations of an isotropic distribution. To compute the expected 
number of random pairs we average over 1, 000 realisations of a 
uniform distribution on the sphere, whose size is given by the num- 
ber of "observed" data points. Note that any small scale clustering 
in angular momentum such as observed in Fig. [T] should manifest 
itself as an excess of pairs with small angular separations. 

Fig. |2] shows the two-point "angular correlation function" 
computed using the present-day angular momentum of subhalos in 
our simulations. Different colours correspond to subhalos accreted 
at different epochs. We see a clear excess of pairs with angular mo- 
mentum orientation separation less than 10° (and up to 30°) com- 
pared to random samples. This implies that the signature of group 
infall is preserved in angular momentum even after many Gyrs of 
evolution. This signal is still discernible even for subhalos accreted 
at z ~ 1. The correlation function calculated with all surviving 
subhalos at present is showed as the (black) dashed line. 

2.2.2 Properties of groups 

We now focus on the characteristics of the groups accreted at the 
various epochs. To identify groups we link pairs of inf ailing halos 
whose angular momentum orientations are separated by a < 10° 
and with relative distances d < 40 kpc at the time of accretion. We 
found that this combination of a and d values results in a robust set 
of groups, maximizing their extent while minimizing the number 
of spurious links. 



3.0 




Figure 2. Two-point "angular correlation function" using the present time 
angular momentum orientation of subhalos accreted in the last 8 Gyrs. The 
excess in the first bins is indicative of the group infall and shows that this 
signal can persist for a very long time. 
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Figure 3. Three examples of the trajectories of groups of subhalos accreted 
at different epochs in the GA3new simulation reference frame. These are 
some of the most abundant groups ever accreted. The colour gradients indi- 
cate the arrow of time, from dark at high redshift to light grey at the present. 
The positions at accretion and present time are highlighted with open circles 
and crosses respectively. There is some evidence that these groups them- 
selves are the result of the mergers of smaller groups, this is especially clear 
for the group in the right-most panel. 



We follow the orbits of the groups identified from redshift z ~ 
4.2 until present time. Fig.[3]shows the trajectories of some of the 
richest groups of subhalos, which were accreted 2.43, 1.65 and 0.84 
Gyrs ago respectively. Each dot represents the position of a subhalo 
colour coded from high-redshift (dark) to the present (light-grey). 
The crosses correspond to the present-day positions while those at 
the time of accretion are shown as open circles. Fig.|3]clearly shows 
that the groups of subhalos follow nearly coherent orbits as early as 
z ~ 4.2, long before the time of accretion. 

The characteristic size of the groups can be measured by com- 
puting the number function of groups, i.e. how many groups have 
a given number of subhalos. Fig. |4] shows the number function of 
groups accreted in the four most recent snapshots: present time, 
0.84, 1.65, and 2.43 Gyrs ago. As can be seen from this Figure, the 
shape is quite similar at all times, and most of the groups have a 
small number of members. 

Fig. |5] shows the differential mass function of the groups in 
Fig. [4] down to our resolution limit (the dashed line, which corre- 
sponds to ~ 5.89 X 10^ Mq). Once again we find very similar 
power-law shapes for the mass functions at different epochs. This 
power-law shape is reminiscent of the differential mass function of 
subhalos in cluster and galaxy-size dark matter halos. The slope of 
the fitted dN/dlogM oc M" relation is n ~ -0.5 ± 0.2. Note 
that this is somewhat shallower than that found for subhalos, where 
n ~ -0.8 ± 0.1 ( |Stoehr et al.p003l [PeLucia et al. 200 4"{ |Gao et| 
|al.|2004bt . This could well be due to insufficient mass resolution: 
the fact that we are not resolving subhalos below 2.9 x 10^ M©, 
implies that many subhalos are accreted in isolation, instead of in 
pairs or in groups. This effect is much stronger at the low mass end 
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Figure 4. Number counts of groups as function of the number of member 
subhalos at four accretion epochs plotted in log-log scale. The error bars are 
Poissonian. 
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Figure 5. Differential mass function of the groups shown in Fig. [4] The 
power-law behaviour of the group mass function is similar to that of the 
subhalo mass function, albeit with a shallower slope. 



of the group mass spectrum. For example, a group with total mass 
~ 10^ Mq can consist of ten subhalos of ~ 10* Mq or two of 
5 X 10* M0. On the other hand a group of ~ lO'^ Mq can only 
be the result of a pair of subhalos of5x lO'' Mq in our simulation. 

Our limited resolution also prevents us from quantifying the 
mass function inside the groups. Nevertheless, and for our largest 
groups we find that these are dominated by a few massive subhalos 
and many small ones. 
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Figure 6. Distribution of angular momentum orientation for the subhalos 
accreted in tlie last 13 snapshots, i.e since z = 1.08. The arrow of time in 
this figure goes from right to left and from bottom to top, i.e. the bottom left 
panel corresponds to ~ 8 ago, while the top left to the present-day. Note 
the small scale structure indicative of group infall as well as the large-scale 
pattern associated to the filamentary structure of the tidal field. This figure 
shows that not just one filament is actively feeding material at any given 
time. 



Figure 7. Evolution of the tidal field around the Milky Way-like halo. The 
grey dots represent 0.05% of the "field" particles within a box of 2/i^^Mpc 
on a side in the simulation reference frame. Like in Fig.|6]the top left panel 
corresponds to the present-time, while the bottom left to 2 ~ 1.08. The 
black circles denote the location of the subhalos accreted at the given epoch. 
The dashed, dash-dot and dotted lines indicate the major, intermediate and 
minor axes of the MW-like halo at each epoch. 



The group infall that we have been detecting in our simulation 
may well be related to the ghostly streams reported by |Lynden-Bell| 
|& Lynden-Bell|(l995^ . The presence of satellites (dwarf galaxies 
and globular clusters) sharing a common orbital plane seems rather 
plausible in the context discussed here. Instead of the disruption 
of a large progenitor ( |Lynden-Bell & Lynden-Belll|1995 < or the 



tidal formation of satellites within gas-rich major mergers < Kroupa| 
|1997 t, we would be witnessing the disruption by the tidal field of 
the Milky Way of a "sub-group"-size object composed by dwarf 
galaxies. The possible implications of this finding are discussed in 
the Conclusions. 



2.2.3 Link to the Environment 

The present-time distribution of angular momentum orientations of 
subhalos reflects both the anisotropy of the accretion pattern and the 
dynamical processes that affect subhalos while orbiting the MW- 
like halo. Fig. [6] shows the orientation of the angular momentum 
of subhalos accreted in the last 13 snapshots, from the present-day 
(top left) to z ~ 1.08 (bottom left). Here the angular momentum is 
calculated using the position and velocity of a subhalo in the simu- 
lation box frame right before it was accreted. Note that only a frac- 
tion of these subhalos will have survived until the present-time. The 
small scale clustering visible in this figure once again highlights the 
group infall. Note the presence of larger-scale patterns lasting over 
several snapshots (in particular in the top row, which corresponds 
to the last 2.4 Gyr). This presumably implies that the infall patterns 
are related with persistent larger scale structures (filaments) in the 
tidal field. 



To understand this in more detail, we proceed to trace the evo- 
lution of the tidal field around the main halo in our simulation. To 
this end we select "field" particles, i.e. those that do not belong to 
the FOF group of the Milky Way-like halo. The projected spatial 
distribution of these particles within a 2/i^^Mpc on a side box is 
shown in grey in Fig. |7] Like in Fig. |6] each panel corresponds to 
a different redshift, starting from z; ~ 1.08 in the bottom left panel 
to the present day in the top left. The distributions of surviving sub- 
halos accreted at the corresponding epoch are overplotted in black. 

Fig. [7] shows that the Milky Way like halo is embedded in a 
larger-scale filamentary pattern. These filaments are comparable in 
extent to the halo itself (as e.g. traced by the accreted subhalos). 
The lumpy nature of the filaments is also clearly visible, showing 
that the infall is not a continuous flow, but is in groups as discussed 
above. 

Note that the global orientation of the tidal fields near the main 
halo has not changed much over the last four snapshots, in agree- 
ment with what is observed in the top row of Fig. |6] Furthermore, 
this large scale pattern is more or less aligned with the major axis of 
the main halo, shown by the dashed line in each panel (as in |Bailin| 
|& Steinmetz|2005| l. 



2.3 On the Great Disk of Milky- Way SateUites 



[Kroupa, Theis & Boily| ( [2005l > and |Metz, Kroupa & Jerjeiil pOOT) ! 
have recently argued that the highly anisotropic distribution of MW 
satellites could not have been drawn from a nearly spherically dis- 
tributed subhalo population. Motivated by their claim and the re- 
sults presented above, we wish to test here under what conditions 
such a configuration is likely in a ACDM simulation like ours. 
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There are many possible ways to define the degree of flatten- 
ing of a distribution. We shaU here concentrate on the following 
two measures: 

(i) The minor- to-major axis ratio c/a derived from the eigen- 
values of the diagonalised inertia tensor defined by the satellites 
positions. 

(ii) The rms of the distances to the best fit plane to the satel- 
lites positions normalized by their median distance from the cen- 
tre: A = Drms/Rmed I Kjoupa, Thcis & Boily|2005 Zentner et 
|al.|2005| i. 

In what follows the positions are defined with respect to the cen- 
troid of the satellites (or subhalos), rather than with respect to the 
centre of the MW(-lik;e) halo. For the first measure (i), we use the 
inertia tensor defined as 



(2) 



Note this inertia tensor differs from that previously used in Eq. 
in which the positions were normalized by their ellipsoidal dis- 
tance. Our preference for this new definition is based on the fact 
that the determination of the ellipsoidal distance is simultaneous 
to the determination of the eigenvalues of the inertia tensor lij. 
This means that an iterative algorithm is used, in which outliers 
are successively discarded, until the desired level of convergence is 
reached (see Sec. |2.2^ . However, when a relatively limited number 
of data points is available (as in the case of the MW satellites) this 
is clearly not desirable. Note as well that the shape of I*j is more 
sensitive to objects at large distances. However, the effect this has 
on the measured c/a can be quantified, as we shall see below. 

In our analysis, we will only consider the eleven "traditional" 
Milky Way satellites. This enables us to make a direct comparison 
to the work of |Kroupa, Theis & Boily| ( f2005[ ). On the other hand, it 
ensures we are not affected by the strong observational bias in the 
sky distribution of the new satellites discovered by SDSS, which re- 
flects the fact that this survey has concentrated on the north galactic 
cap. 

The minor-to-major axis ratio for the set of eleven "tradi- 
tional" MW satellites is c/a ^ 0.18 ± 0.01 where the uncertainty 
is due to errors in the Galactocentric distance (which also includes 
the uncertainty in the distance from the Sun to the Galactic centre, 
Rq = 8.0 ± 0.5 kpc). The plane containing the major and interme- 
diate axes is inclined 72.8 ± 0.7° with respect to the Galactic disk, 
in good agreement with the value found by |Metz, Kroupa & Jerjen| 
( |2007) . 

If we use measure (ii) to quantify the degree of flattening, we 
find that the best-fit plane of the MW satellites has an orientation 
identical to that of the inertia tensor. The rms distance to this plane 
is Drms ~ 18.5 kpc and the plane is offset from the Galactic centre 
by 7.83 kpc. The median distance of the satellites is ~ 80 kpc, 
which implies that A = 0.23 ± 0.01 assuming Gaussian errors of 
all related distances (see also |Metz, Kroupa & Jerjen|2007| >. 

2.3.1 Overall distribution of subhalos 

The subhalos in our GAnew simulations show at least two differ- 
ences in their distributions in comparison to the Galactic dwarfs: 1) 
their spatial distribution is much less anisotropic; 2) their density 
distribution is much shallower. 

The present-day subhalo population of the Milky Way-like 
halo has a minor-to-major axis ratio c/a ~ 0.77, obtained using 
Eq. |2|, compared to c/a = 0.74 for the MW-like halo itself. These 
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Figure 8. Minor-to-major axis ratio of the inertia tensor defined by random 
sets containing a vaiying number of subhalos. The grey dots are for present- 
day subhalos within 300 kpc in GA3new, while the black dots have been 
drawn from an isotropic distribution. The error bars show the 1 — a level 
of each distribution. The c/a obtained by imposing that 11 subhalos follow 
the MW satellites radial distribution is denoted by the triangle, while the 
star symbol corresponds to the c/a for the 11 MW satellites. 



values are fairly consistent with those published by [Zentner et al.| 
(j2005j and |Libeskind et £0^^2005 ). The minor to major axis ratio 
for those subhalos within ~ 300 kpc is c/a ~ 0.86. All these val- 
ues are significantly larger than the c/a ~ 0. 18 found for the Milky 
Way satellites. 

We now wish to test the effect of small number statistics in 
the determination of the principal axes of the inertia tensor. To this 
end we randomly select 10^ samples of A*' subhalos: i) from the 
simulated Milky- Way like halo population; ii) from an isotropic 
distribution with the same radial profile as found in the simulation. 
In Fig. [s] we plot the mean c/a as a function of the number TV of 
subhalos selected. The grey dots correspond to those drawn from 
present-day subhalo population (case i) while the black points are 
drawn from an isotropic sample (case ii). The error bars denote 
the standard deviations of the 10^ realisations. This figure shows 
that only when A'^ ~ 10^ a reliable estimate of the shape of the 
parent distribution can be obtained. The mean c/a tends to become 
smaller, i.e. the distribution appears more flattened as the number of 
subhalos becomes smaller. This is also true when the parent distri- 
bution is completely isotropic. When only 11 subhalos are selected, 
the mean c/a ~ 0.51 ± 0.12 (c/a ~ 0.52 ± 0.12 for the isotropic 
case), quite different from that of the parent sample. Nevertheless, 
this is still significantly larger than observed for the MW satellites. 

We noted previously that the density distribution of the satel- 
lites of the Milky Way is strongly centrally concentrated. The me- 
dian radius of the satellite's distribution is ~ 80 kpc, comparable 
to the half-mass radius of the Milky Way (0.29 r^ir for a concen- 
tration of c = 18 as in |Battaglia et al.|2005^ . On the other hand, the 
median radius for the subhalos (within 300 kpc) in our simulation 
is much larger, Rmed ~ 0.64ri,ir. 
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We now quantify the effect of this highly centrally concen- 
trated distribution on the shape of the inertia tensor. We proceed 
by randomly selecting sets of 11 subhalos from our simulation, but 
now imposing they follow the observed MW satellites spherically 
averaged spatial distribution. The mean minor-to-major axis ratio 
obtained in this way is c/a ~ 0.40 ± 0.10, and is denoted as a 
triangle in Fig. [8] This exercise shows that the disk-like configura- 
tion of the MW satellites is indeed partially driven by their strongly 
centrally concentrated density distribution around the Galaxy (see 
also |Kang et al.|2005||Zentner et al.|2005[ for a similar discussion). 
However, this is still only marginally consistent with the MW satel- 
lites, i.e. the mean c/a of 10^ 11 -random- subhalo realisations is 
~ 2.2(7 away from that observed. 
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2.3.2 Distribution of grouped subhalos 

Inspired by the group infall of subhalos found in Sec. |2.2| we now 
explore how the presence of such groups of subhalos can affect the 
chance of obtaining a flattened configuration. Naively, we would 
expect that if their angular momenta is nearly conserved, the subha- 
los would spread along their orbit and give rise to a planar structure 
as observed. 

Given the large sensitivity of c/a to sample size, we prefer to 
use the A measure in what follows. This measure intuitively ap- 
pears to be more robust since it derives from a simple plane-fit to a 
distribution of points. 

Of the 3, 246 subhalos within 300 kpc from the centre of the 
MW-like halo that have survived until the present-day, 898 subha- 
los fell in as part of a group. 321 different groups have contributed 
to the present-day population of subhalos, of which the earliest two 
were accreted at z — 3.05. From now on, subhalos identified to be 
part of a group are referred to as "grouped", while those which are 
not are termed "field" subhalos. 

We first test how often a disk-like structure is obtained when 
selecting a set of 11 subhalos consisting of a certain number Nsitt 
of subhalos from one group (Ngroup = 1) while the rest N field are 
from the "field". Note that Nsub = 2, 3, . . . , 11, and the condition 
Nsat ~ N field + Nsub = H has to be satisfied. We make lO'' 
realisations of such subsets and compute the fraction that gives rise 
to a configuration as flat as observed, i.e. A 5C 0.23. The result is 
shown in the left panel of Fig. |9] where the disk fraction increases 
from 4.5% to 73% as the number of selected subhalos Nsub in- 
creases from 2 to 11. In comparison, 11-randomly selected subha- 
los (within 300 kpc) gives rise to flattened configurations ~ 2.2% 
of the time. This shows that if the Milky Way satellites fell in to- 
gether, it would not be very surprising that they would be in a planar 
configuration at the present-day. 

It is important to note that when Nsub 6, we select pre- 
dominantly from just two groups accreted at relatively high redshift 
{z = 1.08 and z = 0.84). Other large groups accreted more re- 
cently are still strongly clustered in space, and hence are discarded 
in this exercise since they cannot be considered as a valid repre- 
sentation of the MW satellite population. Furthermore, although 
there is a relatively high chance of obtaining a value of A as low 
as observed, this is driven more by the large median distance of the 
subhalos than by their RMS distance to the best fit plane. 

A second possibility is to consider only "grouped" subhalos. 
That is, we select randomly 11 subhalos from Ngroup different 
groups where Ngroup = 1...5. The panel on the right of Fig. |9] 
shows the fraction of disk-like configurations obtained in this way 
as a function of the number of groups considered. This fraction 
can be as high as ~ 40% when the subhalos come from only two 



Figure 9. Left: Fraction of disk-like structures obtained in 10^ realisations 
consisting Ngub subhalos extracted from one group and 11 — Nsub from 
the "field". Right: Fraction obtained when 11 subhalos are extracted from 
Ngroup different groups. The fraction increases as more subhalos from one 
group are selected, and as the number of groups contributing decreases. The 
likelihood of obtained a highly flattened distribution is always higher than 
when 11 subhalos are randomly selected (dashed lines). 

groups, and of course reaches 73% when they come from just one 
group. Note that the fraction when selecting from 5 different groups 
is still much higher than if one selects 1 1 subhalos randomly. 

The reason for the larger number of disk-like configurations 
when selecting subhalos from several groups as in the right of 
Fig. [9] rather than from just one group and the field, may be un- 
derstood by examining Fig.[TO] This shows the present-day angular 
momentum orientations for "grouped" (left panel) and "field" (right 
panel) subhalos. The two distributions differ clearly in the sense 
that the "grouped" subhalos are generally more clustered (see also 
Sec. \2.2\ , also on larger scales. On the contrary, "field" subhalos 
tend to have their angular momenta more isotropically distributed. 
Therefore, when selecting 11 subhalos purely from groups, the 
chance of picking up subhalos with similar angular momentum ori- 
entations, is higher than when selecting also from the "field". The 
more isotropically distributed orbits of "field" subhalos essentially 
add noise to the highly correlated orbits of subhalos originating in 
just one group. Therefore the disk signal gets smeared out quite 
significantly when more than half of the subhalos are in the field in 
one realisation, as shown on the left of Fig.|9] 

Given the large fraction of flattened configurations found in 
our simulations, we conclude that the spatial distribution of the 
11 Milky Way satellites can be reproduced within ACDM. The 
requirement is that these satellites fell onto the Galactic halo in 
groups. 



3 DISCUSSION AND CONCLUSIONS 

We have analysed a high-resolution cosmological simulation of the 
formation of a Milky Way-like halo, focusing on the properties of 
the satellite population at the present-day. 

We have found that dark-matter subhalos are often accreted in 
groups in our simulations. Roughly 1/3 of the surviving subhalos 
with mass ^ 2.9 x lO'' M0 at the present epoch share this property. 
This is clearly a lower limit since we are not be able to identify 
accompanying halos below our resolution limit (this is particularly 
severe at high-redshift). 

This group infall is apparent as an enhancement in the number 
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Figure 10. Present-time angular momentum orientation of subhalos located 
within 300 kpc from the MW-like host. The left panel shows the subhalos 
from "groups" while the right corresponds to (a randomly selected subset 
from) the "field". 



of subhalos whose angular momentum orientation is similar, par- 
ticularly at the time of infall. This signal is measurable also from 
the present-day angular momentum of subhalos, even for those ac- 
creted 8 Gyr ago. These groups of subhalos share coherent orbits 
which can be traced back well before the accretion epoch. The 
differential group mass function follows a power-law distribution 
dN/d\ogM\ group oc Af" with n ~ -0.5 ± 0.2. This is remi- 
niscent of the differential mass function of subhalos in both galaxy 
and cluster-size halos, albeit with a shallower slope (compared to 
n ~ -0.8 as in e.g. |De Lucia et al.|2004[|Gao et al.|2004b| . 

We have also studied the degree of flattening of the spatial dis- 
tribution of subhalos in our simulation. The mean minor- to-major 
axis ratio c/a of the inertia tensor defined by the positions of 11 
randomly selected subhalos with 300 kpc is c/a ~ 0.51 ± 0.12. 
In comparison the c/a of the 11 "classical" MW satellites is 
0.18 ± 0.01. Imposing the centrally concentrated MW satellite ra- 
dial distribution leads to c/a ~ 0.4 ± 0.1 and therefore somewhat 
alleviates the discrepancy with the observations (see also |Kang et| 
|al.|2005| i. 

We have explored also how this planar configuration may be 
obtained as a result of the infall of satellites in groups. The observed 
correlation in the angular momentum orientation of subhalos nat- 
urally gives rise to disk-like configurations. For example, we find 
that if all subhalos are accreted from just one group, it is almost im- 
possible to avoid a disk-like distribution (~80% probability), while 
for accretion from just two groups, the likelihood of obtaining a 
distribution as planar as observed is 40%. 

These results may explain the origin of the ghostly streams 
proposed by |Lynden-Bell & Lynden-Bell| ([1995). Out of the 
streams originally proposed, only two appear to have survived 
the rigour of time, after modem and accurate measurements of 
proper motions have become available. |Palma, Majewski & John-| 
|ston| ( |2002") confirmed the LMC-SMC-UMi-Draco stream forms a 
kinematic group whose angular momentum separation is < 18.5°. 
More recently [Piatek et al.| ( |2005) ruled out with 95% confident 
level UMi as a member using HST proper motions. The latest mea- 
surements of the Fornax proper motion by Piatek et al.| ^2007^ 
has apparently confirmed the Sculptor-Sextans-Fomax stream (al- 
though previous measurements led to conflicting results, see |Piatek| 
[etaT. 200 2 . Di nescu et al.|2004^ . If some of the luminous satellites 
are embedded in dark (sub)halos that fell in together, such coher- 
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ent structures would be a naturally consequence of the hierarchical 
build-up of galaxies. 

In our simulations, such groups remain coherent in angular 
momentum (i.e. they share similar orbital planes giving rise to great 
circle streams) for approximately 8 Gyr. This implies that these 
groups (or satellites) should have been accreted by the Milky Way 
at redshifts 2 ~ 1 or below. 

One of the possible implications of the reality of the ghostly 
streams is that its member galaxies formed and evolved in a sim- 
ilar environment before falling into the Milky Way potential. This 
would have implications on the (earliest) stellar populations of 
these objects, such as for example, sharing a common metallic- 
ity floor (Helmi et al. 20061. On the other hand, this implies that 
there should be groups that have not been able to host any luminous 
satellites. This would hint at a strong dependence on environment 
on the ability of a subhalo to retain gas ( [Scannapieco et al.'2001 
or be shielded from re-ionization by nearby sources ( Mashchenko, 
[Carignan & Bouchard|2"004l|Weinmann et al.|2007l l. 

Recent proper motion measurements of the Large and Small 
Magellanic clouds by |Kallivayalil, van der M arel & A lcock| ( [2006| l, 
as well as the simulations by Be^T^'cEibal ( )2005j l suggest that 
these systems may have become bound to each other only recently. 
This would be fairly plausible in the context of our results. The 
Clouds may well have been part of a recently accreted group (see 
also [Sales et al.|2007| for the link to Leo I) and it may not even 
be necessary for them to ever have been a binary system. This may 
also have implications on the computations of the past trajectories 
of these systems, particularly if both are embedded in a larger com- 
mon dark-matter envelope ( jBesla et al.|2007^ . 

Our analysis shows that the dynamical peculiarities of the 
Milky Way satellites can be understood in the context of the concor- 
dance cosmological model. Their properties must be a consequence 
of both the environment as well as of the hierarchical nature of the 
build up of galactic halos. 
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